The restoration of the endangered Sambucus palmensis after 30 years of conservation actions in the Garajonay National Park: genetic assessment and niche modeling

Study species

Sambucus palmensis Link, (Sambucus nigra subsp. palmensis (Link) R. Bolli) (Sambucaceae), It grows in the laurel forest vegetation zone between 600–1,000 m a.s.l., with a preference for shady and humid places. It can also be found on the margins of agricultural fields being grown by locals due to its historical use for medicinal purposes (Beltrán et al., 1999; Abdala et al., 2014). This tree can reach up to 4–6 m, is a hermaphroditic taxon and the fruits can be dispersed by birds (Marrero et al., 2011). Also, the species extensively propagates by vegetative reproduction, which has facilitated its reintroduction by cuttings. However, the low regeneration capacity of the species and the high mortality rates of new individuals leads one to think that Sambucus palmensis presents reproductive self-incompatibility and inbreeding (Marrero, Bañares & Carqué, 1998).

Sample collection

The fieldwork was performed under the project 255/2011 funded by “Organismo Autónomo de Parques Nacionales” and the sample collection was assisted by the Garajonay National Park staff. Leaf samples from all the naturally occurring individuals (47), present in the sites of Acebiños, Ancón de Candelaria, Liria, Meriga and Presa de Las Rosas, and a significant representation of the reintroduced individuals from the whole distribution in La Gomera (355) were collected during 2012 and 2013. The categories of natural or reintroduced were considered according to the Garajonay National Park records and annotations. The individuals from Liria and Presa de Las Rosas are all natural, with no reintroduced individuals. In total, 402 samples were collected and genotyped for the study. They are distributed in 15 sites according to their geographical distribution or management area for the National Park. The distribution of the analyzed individuals and the sampling sites are represented in Fig. 1. Young leaves for all specimens were collected and stored in plastic bags with silica gel for their conservation.

Microsatellite development and genotyping

In this study, seven SSR markers were used for the genotyping of Sambucus palmensis individuals in La Gomera. Five markers that were developed for Sambucus nigra (Clarke & Tobutt, 2006) and characterized for S. palmensis (Sosa et al., 2010) were tried for our set of samples. However, only three of these markers (EMSn017, EMSn025, EMSn003) yielded enough polymorphism in the La Gomera populations to be considered. Therefore, we developed four new microsatellite markers for S. palmensis to improve the accuracy of the genotypes, coded as Sam_Tet2, Sam_Hex2, Sam_Hex1 and Sam_Tri8 (Supplemental Information 1).

DNA was extracted from silica-gel-dried young leaves using a modified 2 × CTAB protocol (Doyle & Doyle, 1987). DNA was subsequently purified with the commercial kit Gene Elute PCR Clean-up (Sigma Aldrich, St. Louis, MO, USA). After the characterization of the microsatellites, the whole set of samples was extracted with Invisorb DNA Plant HTS 96 KIT INVISORB.

Microsatellite loci were selected from an Illumina paired-end shotgun library developed by the company AllGenetics (University of A Coruña) using three probe mixes of S. palmensis. We initially chose 20 primer pairs from this library of which 16 yielded some product and were labelled. Finally, four primer pairs amplified consistently with more than two alleles and were selected to complete the genotyping for all samples. For the initial testing, PCR was conducted individually with each primer pair in a 25 uL total reaction volume, which contained approximately: 20 ng of DNA 10 pmol of each primer, as well as PCR Master Mix until 25 uL were completed (Reddy-Mix, ABgene, Surrey, UK). Amplifications were performed using the following conditions: 3 min denaturation at 95 °C, 35 cycles of 30 s denaturation at 95 °C, 30 sannealing at 50–62 °C, and 72 °C for 1.5 min; followed by 5 min elongation at 72 °C. Reverse primers were color-labelled at the 5′-end with 6-FAM, PET, NED or VIC.

Once the new markers where characterized, we conducted the subsequent multiplex amplifications for the seven primer pairs using the QIAGEN Multiplex Kit (QIAGEN). PCR were performed in 15 µL reaction volumes: 7.5 µL of Multiplex PCR Master Mix, 1.5 µL primer mix (containing 0.2 µM of each primer in TE), 1.5 µl of Q-solution, 20–40 ng of DNA and dH₂0. Multiplexing was carried out in two primer groups as indicated in Supplemental Information 1. Following the manufacturer’s instructions, PCRs consisted of a Touchdown protocol with the thermal conditions: 15 min at 95 °C, 10 cycles of 30 s at 94 °C, annealing for 90 s at 65 °C with a decrease of 0.5 °C per cycle and 60 s at 72 °C, followed by 20 cycles of 30 s at 94 °C, annealing for 90 s at 55 °C and 60 sat 72 °C, with a final extension of 30 min at 60 °C. All the products from both simple and multiplex PCR were detected on an ABI 3730 Genetic Analyzer and fragments were sized against the LIZ (500–250) size standard (Applied Biosystems, Foster City, CA, USA) and visualized using Genemapper 4.0 (Applied Biosystems, Foster City, CA, USA). We identified allele peak profiles at each locus and assigned a genotype to each individual.

Statistical analysis

To estimate the incidence of clonality and identify the different genotypes for S. palmensis in La Gomera, multilocus matches were identified with GenAlex version 6.5 (Peakall & Smouse, 2012). Unique genotypes (only in one individual) were coded with consecutive numbers, and shared genotypes (in more than one individual) were coded with letters, from A to ZZ. With the same software, basic genetic diversity indices: average of alleles per locus (Na); number of private alleles (Pa), rare alleles (present in four sites or less; Ra), observed (H_o), and unbiased expected (H_e) heterozygosities for each locus and locality were estimated. Measures of allelic (Ar) and private allelic richness (Par) were calculated using HP-RARE 1.0 (Kalinowski, 2005), which uses rarefaction to correct for sampling error. To detect differences in genetic diversity before and after the restoration programs, the indices calculated for each locality were also estimated for the restored and natural individuals separately. Individual heterozygosity was also calculated to provide a list of the most suitable individuals for restoration programs, in order to increase the fitness of the founder individuals and their offspring (Engelhardt, Lloyd & Neel, 2014). On this list, individuals with unique genotypes, a high individual heterozygosity and presence of rare alleles were identified.

Estimation of null alleles for each locality was carried out with MICROCHECKER 2.2.3 (Van Oosterhout, Weetman & Hutchinson, 2006). BOTTLENECK 1.2.02 software was used to identify any recent genetic drift events in the original set of individuals (Cornuet & Luikart, 1996). The two-phase mutation model (TPM), which is a modification of the stepwise mutation model (SMM), was implemented and shows to be a better fit for most microsatellite data sets (Piry, Luikart & Cornuet, 1999). In the TPM model, to optimize the most sensitive values for microsatellites, a proportion of SMM in the TPM = 0.000 and a variance of the geometric distribution for TPM = 0.36 were chosen.

Allele frequency information was analyzed using a nested analysis of molecular variance (AMOVA) (Excoffier, Smouse & Quattro, 1992) with ARLEQUIN software. The analyses were conducted with two different approaches, with the individuals being grouped either by locality or their origin (natural versus restored). In addition, a principal coordinate analysis (PCoA), using the covariance standardized method of pairwise codominant genotypic distances among individuals, was implemented with GENALEX version 6.5 (Peakall & Smouse, 2012).

To estimate the current genetic structure of the populations, all the genotypes were screened using a Bayesian admixture procedure with the software STRUCTURE (Pritchard, Stephens & Donnelly, 2000). The model was assumed to be of population admixture and correlated allele frequencies. 10 independent runs were conducted for each value of K (from 1 to 15) and analysis consisted of a 10⁵ burn-in period replicated and a run length of 10⁶ replicates. The optimal number of clusters was found by the ΔK method (Evanno, Regnaut & Goudet, 2005) visualized with STRUCTURE HARVESTER (Earl & vonHoldt, 2012). Results of 10 replicates of the best fit K were processed using CLUMPP 1.1.2 (Jakobsson & Rosenberg, 2007) to determine the optimal clustering. The STRUCTURE HARVESTER results for the election of the optimal K are presented in Supplemental Information 2.

Species distribution modeling

The explanatory power of SDMs (Species Distribution Models) is often improved when accurate resolutions of spatial predictors are used (Lassueur, Joost & Randin, 2006; Austin & Van Niel, 2011). We therefore developed a set of accurate spatial climate layers at a 50-meter resolution, based on the network of meteorological stations of the whole archipelago (data provided by AEMET, http://www.aemet.es). Stations with 10 or more years of climate records were filtered for the models. For the monthly variables of minimum and maximum temperature and precipitation, a generalized additive model (GAM) with altitude, northness, latitude and longitude as predictor variables was developed. The best model was selected based on the AIC scores and then projected to La Gomera island. To account for spatial biases of the models, residuals of each model in each meteorological station were used to develop an interpolated map of residuals for each variable by ordinary kriging. This interpolated layer was added to the predicted value of the GAM model to obtain the final layers of each monthly variable. The final dataset of monthly variables was used to develop the bioclimatic variables described by Hijmans et al. (2005), using the ‘dismo’ package in R (Hijmans et al., 2015). Two topographic predictors: slope and topographic index (TPI), derived from the digital elevation model (DEM) of the archipelago were also incorporated.

A total of 441 reintroduced presence cells were recorded in total at 50 meter working resolution, 47 were natural presences in La Gomera and nine recorded in Tenerife and La Palma. A total of 51 presences were discarded as they could not be unambiguously assigned to natural or reintroduced populations in Tenerife and La Palma. To select climatic predictor variables for the modeling procedure first a correlation analysis was conducted for the values of the variables in the cells where the species was present with the R package ‘ecospat’ (Broenniman et al., 2014), which returned a value of four predictor variables. Then a PCA was conducted, and the 11 predictors which obtained the highest scores along the first three axes were retained. Finally, a hierarchical partitioning analysis was conducted (Chevan & Sutherland, 1991), with ‘hier.part’ R package (Walsh & MacNally, 2013), to select among those 11 variables, the five showing the highest independent contributions according to ecospat. These were bioclimatic variables 2, 10, 12 and 15 (temperature range, mean temperature of warmest quarter, annual precipitation, precipitation seasonality, respectively) and slope.

Three different niche modeling procedures were conducted with biomod2 package (Thuiller et al., 2009). The first one was calibrated with natural occurrences only, and the second one with both natural and introduced presences. Finally, the third one was calibrated with introduced occurrences only, to evaluate the model with the natural occurrences and obtaining therefore an evaluation of the introduction with independent data. In each of the three modeling procedures, five datasets were developed containing presence points and pseudoabsences were generated randomly only in unaltered habitats of the island. This was done to avoid generating pseudoabsences in potentially suitable areas were the species was not appearing due to habitat destruction instead of ecological constraints. In the third modeling procedure, pseudoabsences were also generated randomly but avoiding the areas where natural populations existed. In each of the three procedures the number of pseudoabsences differed: with natural occurrences, a total of 200 pseudoabsence points was generated, which were weighted in the models to account for the same importance as presences. With both natural and reintroduced presences, and with only introduced presences, the number of pseudoabsences matched the number of presences.

Six algorithms available in biomod2 were used: generalised linear models (GLM), generalised additive models (GAM), both with stepwise selection, boosting regression (GBM), multiadaptive regression splines (MARS), annual neural networks (ANN), and random forest (RF). 10 runs of each presence-pseudoabsence dataset were run for the first two modeling approaches. In each run, 85% of data was randomly selected for calibration and the rest for model evaluation. The third procedure was evaluated using the natural occurrences of the species. Models were evaluated by means of TSS and ROC scores. Any model with a score below 0.8 for any of the two metrics was excluded for the ensemble model. The remaining models were retained to build an ensemble model based on the contribution of each individual model weighted according to the TSS score.

Species distribution modeling

The PCA analysis conducted with all spatial predictors considered accounted for a 96.35% of the variance and showed patent differences between the niche of introduced and natural populations (Supplemental Information 6), corresponding to differences in the precipitation seasonality (bio04) and annual precipitation (Supplemental Information 7). Also, there were visible contributions of three intertwined variables: precipitation of the wettest month, precipitation of the wettest quarter and precipitation of the coldest month (variables 13,16,19); natural populations showed less temperature seasonality whereas they occurred in areas with lower precipitation (Supplemental Information 7).

The evaluations of the two modeling procedures differed significantly between occurrence datasets (GLMM, F_2,317 = 20.73, P ± 3.575 × 10⁻⁹, pseudoabsence dataset and model run as random effects) and models (GLMM, F_5,317 = 8.5499, P ± 3.575 × 10⁻⁹), with the introduced model showing better performance (Supplemental Information 8). In total, 66 models were not included with the natural dataset, 64 with the introduced one (evaluated against the natural occurrences) and 104 when considering all occurrences. Hence, the ensemble modeling was built with 85, 88, and only 51 respectively under the demanding quality threshold set for filtering models. However, the average TSS score for each modeling dataset was above 0.7, which is considered as fair (natural dataset, TSS = 0.79, introduced 0.78, all = 0.79) (Guisan, Thuiller & Zimmermann, 2017). Models in general showed therefore good performance, yet many of them were not included as they showed values around 0.7–0.8 score.

The suitable areas declared by the three modeling procedures projected suitable areas in the northern parts of the island in areas exposed to humid trade winds, but were conspicuously different (Fig. 3): the suitable areas for S. palmensis appeared in the northern-northeastern section of La Gomera when the model was calibrated with natural presences only. Models calibrated with introduced occurrences only showed on the contrary slightly more extended potential distribution for the species encompassing eastern regions. Regarding the comparison between natural calibration and the model calibrated with all ocurrences, the second showed lower suitability but also a lower threshold. Therefore, the whole model showed a higher extension. The potential distribution fell within protected areas only partially, with 1,345 cells inside protected areas out of 5,045 (26.83%). Only 17 cells did fall within the limits of Garajonay National Park.

Effects of restoration activities on genetic diversity and structure

Because of the main propagation method, by the reintroduction of cuttings, we have detected a high presence of clonal specimens. Just like its relative Sambucus nigra, S. palmensis also presents natural vegetative reproduction (Marrero, Bañares & Carqué, 1998; Bañares et al., 2004). Therefore, natural clones can be found in Liria, which presented 17.14% of shared genotypes. But the percentage of shared genotypes in the other populations that have been restored is much higher, due to the extensive reintroduction of clone specimens. On the other hand, natural regeneration has been detected in some of the sites, with the appearance of new unique genotypes and private alleles in the restored sites. Nonetheless, the private alleles that are present in the restored sites could have come from natural individuals which are now dead and not included in this study. Therefore, the translocation of genotypes may be increasing the chances of sexual reproduction, as has already been detected in some of the restored sites (Marrero, Bañares & Carqué, 2015).

The percentages of variation found between populations are low for an outcrossing species (Hamrick & Godt, 1996), and they are also lower than that found by Sosa et al. (2010) among populations within islands (15%). The variation values detected are also lower than those found for the relative Viburnum treleasi, which has low variation among populations and also presents both sexual and clonal reproduction (Moura, Silva & Caujapé-Castells, 2013). Therefore, there is evidence that the admixture of genotypes across populations have favored the gene flow across the habitat range of S. palmensis in La Gomera. In the STRUCTURE results, Liria, which represents the best-conserved natural population in La Gomera, was the only locality that presented less admixture of individuals. The high number of unique genotypes and the presence of private alleles in Liria suggest that these individuals have been rarely used as a genetic source in the restoration programs, although some individuals were used for propagation before the expansion in 2008 (Á Fernández-López, personal comment, 2017). Although we did not find a sign of recent bottleneck events in the natural source, most of the individuals came from Liria, which has naturally increased its population size since the monitoring programs started. Thus, bottleneck events in the other natural sites are difficult to infer, as only a few individuals have remained.

Overall, there was a light increase in the observed heterozygosity in the restored sites. These results, together with the high admixture found in the genetic structure, also explain the artificial gene flow implemented with the restorations. Outcrossing of inbred isolated populations is playing a major role in the genetic rescue of endangered species (Love Stowell, Pinzone & Martin, 2017), but a balance between genetic rescue and “outbreeding depression” must be found in the management of populations, paying attention to the particular needs in each case (Hufford & Mazer, 2003). In addition, the habitat and climatic continuity of the laurel forest in La Gomera, together with the outcrossing system of the species could have favored the gene flow between the past populations of S. palmensis, hindering high population differentiation or local adaptation. Even if the translocations of genotypes have led to outbreeding depression, the advantages of outcrossing can be greater, especially for self-incompatible species, as it increases the chances of finding available mates (Willi et al., 2007; Pickup & Young, 2008). Despite the increased observed heterozygosity, the number of alleles (the average of alleles per locus and allelic richness) are still lower in the restored individuals compared to the natural ones. These results can be taken as positive, because they indicate that the restored sites can still benefit from a greater outcrossing and admixture with the natural individuals within La Gomera.

Insights from SDMs

This is the first modeling approach at such accurate scale developed for La Gomera. The development of three models based on different datasets for calibration allowed to draw two conclusions about the accuracy of the decision on the placement of the reintroductions. The PCA and the projections of the models clearly show differences in the ecological niche of the natural populations and the introduced ones, which survived successfully under other conditions. Natural populations were localized in areas with lower precipitation and allowed to predict only a part of the total potential area. From the discrepancy between the two niches, three main conclusions can then be drawn.

The first one is that the natural populations, which are very restricted, show a truncated niche, as they only reflected part of the ecological conditions of the species and many individual models did not reach the quality threshold to be included. From the conservation perspective, this is a frequent episode for threatened species witnessing the reduction of the realized niche, and a challenge for modeling procedures which normally should assume niche equilibrium This is the first case in the Canaries of assessment of niche filling of a species through practical reintroduction. Indeed, this kind of assessments supported by independent presence records are rather infrequent (Guisan, Thuiller & Zimmermann, 2017). Therefore, the second conclusion is that the selection criteria of enclaves for reintroduction was appropriate, given the physiological performance of individuals, and despite the climatic conditions did not match those of the scarce natural remaining populations. The quality of the models that considered introduced occurrences, assessed either with independent data or not, was significantly higher than considering natural occurrences only, with few models failing. The fairly higher number of introduced occurrences on the other hand can explain why the potential distribution predicted by the model calibrated with introduced occurrences and absences within La Gomera is higher than any other model. Finally, models served as a tool to identify further areas for protection. At present, protected areas other than Garajonay National only cover a quarter of the potential distribution.

Recommendations for conservation actions

One of the main purposes of restoration ecology is to simulate the characteristics of the natural populations (Pavlik, 1996) and the restored sites do not show diminished levels of genetic diversity compared to the original populations, despite the high number of clonal specimens. Also, SDMs support that reintroduced specimens match properly missing ecological conditions in the remaining populations. However, to improve the sexual regeneration in future reintroductions, further studies of the reproductive biology of S. palmensis are encouraged. The detection of the possible causes of self-incompatibility would help to increase the level of available mates and therefore gene flow and offspring. Also, more demographic studies such as that carried out in Meseta Vallehermoso (Marrero, Bañares & Carqué, 2015), will help to monitor the fitness and survival of restored sites over time. The combination of demographic and genetic studies is vital to ensure the recovery of endangered species (Oostermeijer, Luijten & Den Nijs, 2003).

As an urgent measure to maintain the genetic diversity of S. palmensis in La Gomera, we have already provided a list to the Garajonay National Park managers with the best candidates for a conservation genetic program (Supplemental Information 5). We also believe that it is important to avoid the genotypes detected in this study that have been extensively used in some sites. As indicated in Vergeer, Sonderen & Ouborg (2004), we would also suggest increasing the number of unrelated seed producer individuals to create sustainable and viable populations, which would also avoid inbreeding processes. Since propagation by cuttings is the most viable way of reintroducing new individuals, the consideration of all individuals with unique genotypes for future reintroductions is also a conservation measure to be taken into account. Although it is possible that individuals from Tenerife were introduced in the past (Sosa et al., 2010), we consider that there is enough genetic variability in La Gomera to continue the restoration programs using the genotypes that are currently present on the island. Moreover, the suitability model based on the whole dataset of occurrences is a valid tool to identify other suitable areas for further reintroductions or translocations. The fact that the current S. palmensis’ populations in La Gomera are within a National Park figure, clearly has benefit the success of the restoration programs, in respect to the populations in other islands. Although we have found suitable areas out of Garajonay that could be useful in future reintroduction attempts. Therefore, it is important to take into consideration endemic rarity in the design of protected areas (Irl et al., 2017).

It has been detected that good germination rates and seedling establishment highly depends on years of good precipitation (Marrero, Bañares & Carqué, 2015). But the long survival of individuals is determined by the clearance of the forest, reason why the Liria population expanded after the fire in 2008. Although introduced herbivores are a major threat for the island endemic flora (Reaser et al., 2007; Garzón-Machado et al., 2010; Irl et al., 2014), the laurel forest is not inhabited by rabbits, and the Garajonay National Park tries to control the effects of feral goats and sheep. However, predation of the fruits and leaves by rats has been observed (Á Fernández, personal comment, 2017). In this sense, we consider of importance to favor the establishment of new seedlings by manually clearing the canopy forest in the surroundings of S. palmensis individuals and a thorough control of invasive herbivores.

On a long-term basis, this case study will provide a great deal of information regarding the consequences of restoration actions in self-incompatible clonal species. Moreover, these results can serve as a guideline for the restoration programs in Tenerife, La Palma and Gran Canaria, whose island governments are also restoring S. palmensis populations. For example, in Gran Canaria, only two naturally occurring individuals and some cultivars were found prior to the reintroductions. Therefore, a better knowledge of the genetic background of the restored individuals and the climatic suitability of the species would increase the success and long survival of the populations.